##
## the year to work with: 
cy<-2005


## step 1: read in pre-prepared data: 

dd <- c(paste("C:/xuncao/mcmc_net_polity/", as.character(cy), "try 2", sep="")) 
setwd(dd)

dat<-dget(paste("data_wd_nomissing_polity", as.character(cy), sep=""))

## for Xd: we have, in exact order, distance between capital cities (-), common language (+), and correlation of growth rate (-) --- expected sign in the (). 
## I am going to replace common language with UTV later on; so common language won't enter the analysis. 

Xs<-na.omit(dat[[3]]) # sender effects:  polity,  corruption (notice that higher values for less corrupted governments); CapOpen; GDP (logged); GDPPC (logged)
Xs[, 4]<-log(Xs[, 4])
Xs[, 5]<-log(Xs[, 5])

Xr<-Xs
# Xr<-na.omit(dat[[3]]) # receiver effects: corruption; capital market opennes; GDP;  GDP per caipta;

Yd<-dat[[1]][rownames(Xs), rownames(Xs)]
Yd[is.na(Yd)]<-0 # be consistent with Xd
                 # I think gbme can do the missing value imputation for Yd; but it didn't work this time.
Yd[Yd<0]<-0      # this is weired from the data, there is a negative value there!!!

diag(Yd)<-NA
Yd<-log(Yd+1)

Xd<-dat[[2]][rownames(Xs), rownames(Xs), ] 




## step 2: use the latent space estimated from 2001: this is the starting year
## the following years actually only estimate the changes in the latent space.

OUT=read.table("C:/xuncao/mcmc_net_polity/2001try 2/OUT",header=T) 
U=read.table("C:/xuncao/mcmc_net_polity/2001try 2/U")
V=read.table("C:/xuncao/mcmc_net_polity/2001try 2/V")
 
nss=dim(OUT)[1]
n=dim(U)[1]/nss
k<-3
nburn<-100
aa<-colnames(Yd)

PU<-array(dim=c(n,k,nss))
for(i in 1:nss) { PU[,,i]<-as.matrix(U[ ((i-1)*n+1):(i*n) ,])  }
PU<-PU[,,-(1:nburn)];rownames(PU)=aa;     #drop burn in
PV<-array(dim=c(n,k,nss))
for(i in 1:nss) { PV[,,i]<-as.matrix(V[ ((i-1)*n+1):(i*n) ,])  }

# drop burn in ##
PV<-PV[,,-(1:nburn)]; rownames(PV)=aa; 
OUT=OUT[-(1:nburn),]

 
#find posterior mean of Z %*% t(Z)
UTV<-matrix(0,n,n)
for(i in 1:dim(PU)[3] ) { ## Do not change k here
    UTV<-UTV+PU[,,i]%*%t(PV[,,i]) 
}
UTV<-UTV/dim(PU)[3]  ## Do not change to k here
tmp<-svd(UTV)
 
u2=tmp$u[,1:k];rownames(u2)=aa;
v2=tmp$v[,1:k];rownames(v2)=aa;
# u2[aa1,]=u1;v2[aa1,]=v1;d2=tmp$d[1:k];


##
# no need to do the following if we have consistent set of countries over the years. 

#Xd=array(dim=c(length(Xs),length(Xs),2));rownames(Xd)=aa
#rownames(Xs)=rownames(Y)
#Xr<-Xs

#    Xd[,,1]=u2%*%diag(d2)%*%t(v2)
#    Xd[,,2]=xd




## step 3: we also need the change in latent space estimated from 2002-to the previous year
## I used U.d and V.d for change; this is cumulative ...

for (i in 2002:(cy-1)){ 
OUT=read.table(paste("C:/xuncao/mcmc_net_polity/", as.character(i), "try 2/OUT",sep=""), header=T) 
U.d=read.table(paste("C:/xuncao/mcmc_net_polity/", as.character(i), "try 2/U",sep=""))
V.d=read.table(paste("C:/xuncao/mcmc_net_polity/", as.character(i), "try 2/V",sep=""))
 
nss=dim(OUT)[1]
n=dim(U.d)[1]/nss
k<-3
nburn<-100
aa<-colnames(Yd)

PU.d<-array(dim=c(n,k,nss))
for(i in 1:nss) { PU.d[,,i]<-as.matrix(U.d[ ((i-1)*n+1):(i*n) ,])  }
PU.d<-PU.d[,,-(1:nburn)];rownames(PU.d)=aa;     

PV.d<-array(dim=c(n,k,nss))
for(i in 1:nss) { PV.d[,,i]<-as.matrix(V.d[ ((i-1)*n+1):(i*n) ,])  }
PV.d<-PV.d[,,-(1:nburn)]; rownames(PV.d)=aa; 

OUT=OUT[-(1:nburn),]


#find posterior mean of Z %*% t(Z)
UTV.d<-matrix(0,n,n)
for(i in 1:dim(PU.d)[3] ) { ## Do not change k here
    UTV.d<-UTV.d+PU.d[,,i]%*%t(PV.d[,,i]) 
}
UTV.d<-UTV.d/dim(PU.d)[3]  ## Do not change to k here
tmp<-svd(UTV.d)
 
u2.d=tmp$u[,1:k];rownames(u2.d)=aa;
v2.d=tmp$v[,1:k];rownames(v2.d)=aa;

UTV<-UTV+UTV.d 
}


## step 4: 
## we need to plug the sum of latent space estimated from year 2001 plus "the change" in latent space estimated from 2002.
Xd[,,2]<-UTV   


## step 5: initiate gbme:
source("gbme.asym.r")
n<-dim(Yd)[1]
Y<-Yd

# call the main estimation package, gbme
gbme(Y=Yd,Xd=Xd,Xs=Xs,Xr=Xr,k=3,NS=100000,awrite=T,bwrite=T, oround=8) 
